extract_coefficients <- function(models, unique_variables) {
  results <- data.frame(variable = rep(unique_variables, each = length(models)),
                        coefficients = numeric(length(models) * length(unique_variables)),
                        standard_errors = numeric(length(models) * length(unique_variables)))
  
  for (i in seq_along(models)) {
    model_summary <- summary(models[[i]])
    for (j in seq_along(unique_variables)) {
      variable_name <- unique_variables[j]
      if (variable_name %in% rownames(model_summary$coefficients)) {
        coef_index <- match(variable_name, rownames(model_summary$coefficients))
        results[(i - 1) * length(unique_variables) + j, ] <- c(variable_name,
                                                               model_summary$coefficients[coef_index, "Estimate"],
                                                               model_summary$coefficients[coef_index, "Std. Error"])
      } else {
        results[(i - 1) * length(unique_variables) + j, ] <- c(variable_name, NA, NA)
      }
    }
  }
  return(results)
}


plot_results <- function(results, x_labels, y_label, colors, file_path, color_column) {
  # Determine the color aesthetic based on the color_column
  if (color_column == "party") {
    color_aesthetic <- aes(x = model, y = coefficients, color = party)
    legend_position <- "bottom"
    color_lab <- 'Party'
  } else {
    color_aesthetic <- aes(x = model, y = coefficients, color = column_set)
    legend_position <- "none"
    color_lab <- "Treatment"
  }
  
  fig <- ggplot(results) +
    color_aesthetic +
    geom_point(position = position_dodge(width = 0.3), show.legend = TRUE) +
    geom_errorbar(aes(ymin = coefficients - 1.96 * standard_errors, ymax = coefficients + 1.96 * standard_errors),
                  position = position_dodge(width = 0.3), width = 0.2) +
    geom_hline(yintercept = 0, color = "black", linetype = "solid") +  
    labs(x = "", y = y_label, color = color_lab) +
    scale_color_manual(values = colors) +   
    scale_x_discrete(labels = x_labels) +
    theme_minimal() +
    theme(axis.text.x = element_text(hjust = 0.5), legend.position = legend_position, plot.caption = element_text(color = "grey27", size = 8)) +
    facet_grid(. ~ column_set, scales = "free_x", space = 'free', switch = "both")
  print(fig)
  ggsave(file_path, plot = fig, width = 7, height = 5)
}


create_results_table <- function(t_test_message_main, t_test_c_donate_1, control_group) {
  tidy_message_main <- tidy(t_test_message_main)
  tidy_c_donate_1 <- tidy(t_test_c_donate_1)
  
  results_table <- bind_rows(
    tidy_message_main %>% mutate(Outcome = "message_main"),
    tidy_c_donate_1 %>% mutate(Outcome = "c_donate_1")
  ) %>%
    dplyr::select(Outcome, estimate, statistic, p.value, conf.low, conf.high)
  
  results_table$mean.black <- c(mean(control_group[control_group$black == 1, ]$message_main, na.rm = TRUE),
                                mean(control_group[control_group$black == 1, ]$c_donate_1, na.rm = TRUE))
  
  results_table$mean.white <- c(mean(control_group[control_group$white == 1, ]$message_main, na.rm = TRUE),
                                mean(control_group[control_group$white == 1, ]$c_donate_1, na.rm = TRUE))
  
  return(results_table)
}

# Function to perform randomization inference
randomization_inference <- function(data, outcome_var, frame_var, level1, level2, n_perm = 10000) {
  
  # Subset data for the specified levels
  data_subset <- subset(data, get(frame_var) %in% c(level1, level2))
  
  # Fit the original model
  formula <- as.formula(paste(outcome_var, "~ factor(", frame_var, ")"))
  main_results <- lm_robust(formula, data = data_subset)
  
  # Determine the coefficient name based on levels
  coef_name <- paste0("factor(", frame_var, ")", level2)
  
  # Extract the observed coefficient
  obs_coef <- coef(main_results)[coef_name]
  
  # Initialize vector to store permutation results
  perm_coefs <- numeric(n_perm)
  
  for (i in 1:n_perm) {
    # Shuffle the treatment labels
    permuted_frame <- sample(data_subset[[frame_var]])
    # Refit the model with permuted treatment
    perm_results <- lm_robust(as.formula(paste(outcome_var, "~ factor(permuted_frame)")), data = data_subset)
    coef_name2 <- paste0("factor(permuted_frame)", level2)
    # Extract the coefficient for the permuted treatment
    perm_coefs[i] <- coef(perm_results)[coef_name2]
  }
  
  # Calculate the p-value
  p_value <- mean(abs(perm_coefs) >= abs(obs_coef))
  
  # Output results
  list(observed_coef = obs_coef, p_value = p_value)
}

